home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_4 / a4_1a.m next >
Encoding:
Text File  |  1994-06-05  |  4.3 KB  |  177 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 4.1 (Evaluation of a Taylor Series).
  9. % Section    4.1,    Taylor Series and Calculation of Functions, Page 203
  10. echo on; clc; format long; hold off; clear
  11.  
  12. % This program investigatges Taylor polynomial approximations.
  13.  
  14. %      Pn(x) = c(1) + c(2)x + c(2)x^2 + ... + c(n+1)x^n
  15.  
  16. % The coefficient list for the Taylor polynomial of degree
  17.  
  18. % n = 5  for  f(x) = sin(x)  is stored in the vector  C
  19.  
  20. % The function  sin(x)  is stored in the string variable  fun
  21.  
  22. C = [1/120,0,-1/6,0,1,0];
  23.  
  24. fun = 'sin(x)';
  25.  
  26. pause    % Press any key to graph f(x) and Pn(x).
  27.  
  28. clc; clg
  29. a = -0.1;
  30. b =  3.2;
  31. c = -0.1;
  32. d =  1.05;
  33. n = length(C)-1;
  34. X = 0:0.01:pi;
  35. x = X;
  36. Y = eval(fun);
  37. P = polyval(C,X);
  38. axis([a b c d]);...
  39. plot(X,Y,'-g',X,P,'-r');...
  40. hold on;...
  41. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  42. xlabel('x');...
  43. ylabel('y');...
  44. Mx1 = ['Comparison of ',fun,' and P'];...
  45. Mx2 = [Mx1,num2str(n),'(x) over ['];...
  46. Mx3 = [Mx2,num2str(a),',',num2str(b),'].'];...
  47. title(Mx3);...
  48. hold off;...
  49. axis;...
  50. shg; pause    % Press any key to continue.
  51.  
  52. Mx1 = 'The function is  f(x) = ';
  53. Mx2 = 'Pn(x) = c(1)x^n + c(2)x^(n-1) + ... + c(n)x + c(n+1)';
  54. Mx3 = 'The degree is  n = ';
  55. Mx4 = ', and the coefficients list  c  is:';
  56. clc,echo off,format short,diary output,...
  57. disp(''),disp([Mx1,fun]),...
  58. disp(Mx2),disp([Mx3,num2str(n),Mx4]),...
  59. for i=1:5:n+1, disp(C([i:min(i+4,n+1)])); end,...
  60. diary off,echo on,format long
  61.  
  62. pause    % Press any key to graph f(x) - Pn(x).
  63.  
  64. clc; clg;
  65. a = -0.1;
  66. b = 1.0;
  67. c = -0.0002;
  68. d = 0.00001;
  69. axis([a b c d]);...
  70. plot(X,Y-P,'-r');...
  71. hold on;...
  72. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  73. xlabel('x');...
  74. ylabel('y');...
  75. Mx1 = ['The error: ',fun,' - P'];...
  76. Mx2 = [Mx1,num2str(n),'(x) over ['];...
  77. Mx3 = [Mx2,num2str(a),',',num2str(b),'].'];...
  78. title(Mx3);...
  79. hold off;...
  80. axis;...
  81. shg; pause    % Press any key for a list of numerical computations.
  82.  
  83. clc;
  84. X = 0:0.05:1;
  85. x = X;
  86. Y = eval(fun);
  87. P = polyval(C,X);
  88. points = [X;Y;P;Y-P];
  89. Mx1=['Taylor polynomial approximation of f(x) = ',fun];
  90. Mx2='     x(k)               f(x(k))            Pn(x(k))           error';
  91. clc,echo off,diary output,...
  92. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(points'),...
  93. diary off, echo on
  94.  
  95. pause    % Press any key for another Taylor polynomial.
  96.  
  97. clc;
  98. % The coefficient list for the Taylor polynomial of degree
  99.  
  100. % n = 9  for  f(x) = exp(x)  is stored in the vector  C
  101.  
  102. % The function  exp(x)  is stored in the string variable  fun
  103.  
  104. C = [1/362880,1/40320,1/5040,1/720,1/120,1/24,1/6,1/2,1,1];
  105.  
  106. fun = 'exp(x)';
  107.  
  108. X = -8:0.05:8;
  109. x = X;
  110. Y = eval(fun);
  111. P = polyval(C,X);
  112.  
  113. pause    % Press any key to graph f(x) and Pn(x).
  114.  
  115. clc; clg;
  116. a = -8;
  117. b =  8;
  118. c = -200;
  119. d = 2000;
  120. n = length(C)-1;
  121. axis([a b c d]);...
  122. plot(X,Y,'-g',X,P,'-r');...
  123. hold on;...
  124. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  125. xlabel('x');...
  126. ylabel('y');...
  127. Mx1 = ['Comparison of ',fun,' and P'];...
  128. Mx2 = [Mx1,num2str(n),'(x) over ['];...
  129. Mx3 = [Mx2,num2str(a),',',num2str(b),'].'];...
  130. title(Mx3);...
  131. hold off;...
  132. axis;...
  133. shg; pause    % Press any key to continue.
  134.  
  135. Mx1 = 'The function is  f(x) = ';
  136. Mx2 = 'Pn(x) = c(1)x^n + c(2)x^(n-1) + ... + c(n)x + c(n+1)';
  137. Mx3 = 'The degree is  n = ';
  138. Mx4 = ', and the coefficients list  c  is:';
  139. clc,echo off,format short,diary output,...
  140. disp(''),disp([Mx1,fun]),...
  141. disp(Mx2),disp([Mx3,num2str(n),Mx4]),...
  142. for i=1:5:n+1, disp(C([i:min(i+4,n+1)])); end,...
  143. diary off, echo on,format long
  144.  
  145. pause    % Press any key to graph f(x) - Pn(x).
  146.  
  147. clc; clg;
  148. a = -1.0;
  149. b =  1.0;
  150. c = -0.00000001;
  151. d =  0.0000003;
  152. axis([a b c d]);...
  153. plot(X,Y-P,'-r');...
  154. hold on;...
  155. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  156. xlabel('x');...
  157. ylabel('y');...
  158. Mx1 = ['The error: ',fun,' - P'];...
  159. Mx2 = [Mx1,num2str(n),'(x) over ['];...
  160. Mx3 = [Mx2,num2str(a),',',num2str(b),'].'];...
  161. title(Mx3);...
  162. hold off;...
  163. axis;...
  164. shg; pause    % Press any key for a list of numerical computations.
  165.  
  166. clc;
  167. X = -1:0.1:1;
  168. x = X;
  169. Y = eval(fun);
  170. P = polyval(C,X);
  171. points = [X;Y;P;Y-P];
  172. Mx1=['Taylor polynomial approximation of f(x) = ',fun];
  173. Mx2='     x(k)               f(x(k))            Pn(x(k))           error';
  174. clc,echo off,diary output,...
  175. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(points'),...
  176. diary off, echo on
  177.